#VERSION INFORMATION
#RESULTS COMPUTED WITH R 3.0.0.
#VERSION 1.01 OF THE 'TSA' LIBRARY
#VERSION 0.10-31 OF THE 'tseries' LIBRARY
#VERSION 1.5-9.1 OF THE 'locfit' LIBRARY
#VERSION 1.7-22. OF THE 'mgcv' LIBRARY
#VERSION 0.8-53 OF THE 'foreign' LIBRARY
#VERSION 1.4-0 OF THE 'plm' LIBRARY
#VERSION 0.2-5 OF THE 'XLConnect' LIBRARY
#In the unlikely event that results are not exactly reproduced with the code in this file, you may try the file versions listed above.

#libraries
rm(list=ls())
library(TSA)
library(foreign)
library(lattice)
library(plm)
library(XLConnect)
library(reshape)

#options
options(digits=10,scipen=16)

###set to YOUR working directory
setwd('/Volumes/MONOGAN/CONSULTING/japanEnergy/civilSociety')

#load data 
POLICYII <- read.dta('japanCityElectric.dta')
summary(POLICYII)
long.series<-readWorksheetFromFile("2001-2012 data_rev.xlsx", sheet=1, header=TRUE)
summary(long.series)

#data cleaning
POLICYII$YearDum<-as.factor(POLICYII$year)
long.stacked.series<-melt.data.frame(data=long.series, measure.vars=c("Sapporo","Sendai","Saitama.City.","Chiba.","Tokyo","Kawasaki","Yokohama.","Niigata","Shizuoka","Hamamatsu","Nagoya.","Kyoto","Osaka","Sakai","Kobe","Hiroshima","Kitakyushu","Fukuoka"),variable_name="city", id=c("Col1"))

###PERCENTAGE DROP IN ELECTRICITY###
POLICYII.b<-subset(POLICYII,city!="Shizuoka ")
(sum(POLICYII.b$dependent[POLICYII.b$year==2012])-sum(POLICYII.b$dependent[POLICYII.b$year==2010]))
(sum(POLICYII.b$dependent[POLICYII.b$year==2012])-sum(POLICYII.b$dependent[POLICYII.b$year==2010]))/sum(POLICYII.b$dependent[POLICYII.b$year==2010])

###VIEW DATA###
#Line graph of electricity consumption by city, 2001-2012
trellis.device("png",color=FALSE,file="cityElectricityLONG.png",height=600,width=600,pointsize=16)
xyplot(log(value)~Col1, data=long.stacked.series, type='l', groups=city, xlab="Year", ylab="Electicity Consumption in GWh (Logarithmic Scale)", scales=list(y=list(at=c(log(seq(from=1000,to=35000,by=1000))),labels=c("","2000","3000","4000","5000",rep("",4),"10000",rep("",4),"15000",rep("",4),"20000",rep("",4),"25000",rep("",4),"30000",rep("",4),"35000"))))
dev.off()

trellis.device("png",color=FALSE,file="cityElectricityNoTokyoLONG.png",height=600,width=600,pointsize=16)
xyplot(log(value)~Col1, data=long.stacked.series, subset=city!="Tokyo", type='l', groups=city, xlab="Year", ylab="Electicity Consumption in GWh (Logarithmic Scale)", scales=list(y=list(at=c(log(seq(from=1000,to=35000,by=1000))),labels=c("","2000","3000","4000","5000",rep("",4),"10000",rep("",4),"15000",rep("",4),"20000",rep("",4),"25000",rep("",4),"30000",rep("",4),"35000"))))
dev.off()

#Line graph of electricity consumption by city, 2007-2012
trellis.device("png",color=FALSE,file="cityElectricity.png",height=600,width=600,pointsize=16)
xyplot(log(dependent)~year, data=POLICYII, type='l', groups=city, xlab="Year", ylab="Electicity Consumption in GWh (Logarithmic Scale)", scales=list(y=list(at=c(log(seq(from=1000,to=35000,by=1000))),labels=c("","2000","3000","4000","5000",rep("",4),"10000",rep("",4),"15000",rep("",4),"20000",rep("",4),"25000",rep("",4),"30000",rep("",4),"35000"))))
dev.off()

trellis.device("png",color=FALSE,file="cityElectricityNoTokyo.png",height=600,width=600,pointsize=16)
xyplot(log(dependent)~year, data=POLICYII, subset=city!="Tokyo", type='l', groups=city, xlab="Year", ylab="Electicity Consumption in GWh (Logarithmic Scale)", scales=list(y=list(at=c(log(seq(from=1000,to=35000,by=1000))),labels=c("","2000","3000","4000","5000",rep("",4),"10000",rep("",4),"15000",rep("",4),"20000",rep("",4),"25000",rep("",4),"30000",rep("",4),"35000"))))
dev.off()

###RANDOM EFFECTS MODELS OF THE DISASTER'S EFFECT ON AVERAGE###
re1=plm(dependent~psdv,data=POLICYII,index=c("city","year"),model="random");summary(re1)

re2=plm(dependent~psdv*npo,data=POLICYII,index=c("city","year"),model="random");summary(re2)

re3=plm(dependent~psdv*enpo,data=POLICYII,index=c("city","year"),model="random");summary(re3)

re4=plm(dependent~psdv*gnpo,data=POLICYII,index=c("city","year"),model="random");summary(re4)

###RANDOM EFFECTS MODELS OF THE DISASTER'S EFFECT CONTINGENT ON POLICY###
re5=plm(dependent~ttdv+kcdv+nopolicy,data=POLICYII,index=c("city","year"),model="random");summary(re5)

re6=plm(dependent~ttdv*npo+kcdv*npo+nopolicy*npo,data=POLICYII,index=c("city","year"),model="random");summary(re6)

re7=plm(dependent~ttdv*enpo+kcdv*enpo+nopolicy*enpo,data=POLICYII,index=c("city","year"),model="random");summary(re7)

re8=plm(dependent~ttdv*gnpo+kcdv*gnpo+nopolicy*gnpo,data=POLICYII,index=c("city","year"),model="random");summary(re8)

###POST HOC T-RATIOS###
#Model 2
#Minimium
effect<-re2$coefficients[2]+(re2$coefficients[4]*min(POLICYII$npo)); effect
var<-re2$vcov[2,2]+2*min(POLICYII$npo)*re2$vcov[2,4]+(min(POLICYII$npo)^2)*re2$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Mean
effect<-re2$coefficients[2]+(re2$coefficients[4]*mean(POLICYII$npo)); effect
var<-re2$vcov[2,2]+2*mean(POLICYII$npo)*re2$vcov[2,4]+(mean(POLICYII$npo)^2)*re2$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-re2$coefficients[2]+(re2$coefficients[4]*max(POLICYII$npo)); effect
var<-re2$vcov[2,2]+2*max(POLICYII$npo)*re2$vcov[2,4]+(max(POLICYII$npo)^2)*re2$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 3
#Minimium
effect<-re3$coefficients[2]+(re3$coefficients[4]*min(POLICYII$enpo)); effect
var<-re3$vcov[2,2]+2*min(POLICYII$enpo)*re3$vcov[2,4]+(min(POLICYII$enpo)^2)*re3$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Mean
effect<-re3$coefficients[2]+(re3$coefficients[4]*mean(POLICYII$enpo)); effect
var<-re3$vcov[2,2]+2*mean(POLICYII$enpo)*re3$vcov[2,4]+(mean(POLICYII$enpo)^2)*re3$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-re3$coefficients[2]+(re3$coefficients[4]*max(POLICYII$enpo)); effect
var<-re3$vcov[2,2]+2*max(POLICYII$enpo)*re3$vcov[2,4]+(max(POLICYII$enpo)^2)*re3$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 4
#Minimium
effect<-re4$coefficients[2]+(re4$coefficients[4]*min(POLICYII$gnpo)); effect
var<-re4$vcov[2,2]+2*min(POLICYII$gnpo)*re4$vcov[2,4]+(min(POLICYII$gnpo)^2)*re4$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Mean
effect<-re4$coefficients[2]+(re4$coefficients[4]*mean(POLICYII$gnpo)); effect
var<-re4$vcov[2,2]+2*mean(POLICYII$gnpo)*re4$vcov[2,4]+(mean(POLICYII$gnpo)^2)*re4$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-re4$coefficients[2]+(re4$coefficients[4]*max(POLICYII$gnpo)); effect
var<-re4$vcov[2,2]+2*max(POLICYII$gnpo)*re4$vcov[2,4]+(max(POLICYII$gnpo)^2)*re4$vcov[4,4]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 6
#Minimium
effect<-re6$coefficients[2]+(re6$coefficients[6]*min(POLICYII$npo)); effect
var<-re6$vcov[2,2]+2*min(POLICYII$npo)*re6$vcov[2,6]+(min(POLICYII$npo)^2)*re6$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Mean
effect<-re6$coefficients[2]+(re6$coefficients[6]*mean(POLICYII$npo)); effect
var<-re6$vcov[2,2]+2*mean(POLICYII$npo)*re6$vcov[2,6]+(mean(POLICYII$npo)^2)*re6$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-re6$coefficients[2]+(re6$coefficients[6]*max(POLICYII$npo)); effect
var<-re6$vcov[2,2]+2*max(POLICYII$npo)*re6$vcov[2,6]+(max(POLICYII$npo)^2)*re6$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 7
#Minimium
effect<-re7$coefficients[2]+(re7$coefficients[6]*min(POLICYII$enpo)); effect
var<-re7$vcov[2,2]+2*min(POLICYII$enpo)*re7$vcov[2,6]+(min(POLICYII$enpo)^2)*re7$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Mean
effect<-re7$coefficients[2]+(re7$coefficients[6]*mean(POLICYII$enpo)); effect
var<-re7$vcov[2,2]+2*mean(POLICYII$enpo)*re7$vcov[2,6]+(mean(POLICYII$enpo)^2)*re7$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-re7$coefficients[2]+(re7$coefficients[6]*max(POLICYII$enpo)); effect
var<-re7$vcov[2,2]+2*max(POLICYII$enpo)*re7$vcov[2,6]+(max(POLICYII$enpo)^2)*re7$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 8
#Minimium
effect<-re8$coefficients[2]+(re8$coefficients[6]*min(POLICYII$gnpo)); effect
var<-re8$vcov[2,2]+2*min(POLICYII$gnpo)*re8$vcov[2,6]+(min(POLICYII$gnpo)^2)*re8$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Mean
effect<-re8$coefficients[2]+(re8$coefficients[6]*mean(POLICYII$gnpo)); effect
var<-re8$vcov[2,2]+2*mean(POLICYII$gnpo)*re8$vcov[2,6]+(mean(POLICYII$gnpo)^2)*re8$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-re8$coefficients[2]+(re8$coefficients[6]*max(POLICYII$gnpo)); effect
var<-re8$vcov[2,2]+2*max(POLICYII$gnpo)*re8$vcov[2,6]+(max(POLICYII$gnpo)^2)*re8$vcov[6,6]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

###CREATE STANDARDIZED VARIABLES###
POLICYII.std<-POLICYII
POLICYII.std$dependent<-(POLICYII.std$dependent-mean(POLICYII.std$dependent,na.rm=T))/sd(POLICYII.std$dependent,na.rm=T)
POLICYII.std$npo<-(POLICYII.std$npo-mean(POLICYII.std$npo,na.rm=T))/sd(POLICYII.std$npo,na.rm=T)
POLICYII.std$enpo<-(POLICYII.std$enpo-mean(POLICYII.std$enpo,na.rm=T))/sd(POLICYII.std$enpo,na.rm=T)
POLICYII.std$gnpo<-(POLICYII.std$gnpo-mean(POLICYII.std$gnpo,na.rm=T))/sd(POLICYII.std$gnpo,na.rm=T)
POLICYII.std$psdv<-(POLICYII.std$psdv-mean(POLICYII.std$psdv,na.rm=T))/sd(POLICYII.std$psdv,na.rm=T)
POLICYII.std$ttdv<-(POLICYII.std$ttdv-mean(POLICYII.std$ttdv,na.rm=T))/sd(POLICYII.std$ttdv,na.rm=T)
POLICYII.std$kcdv<-(POLICYII.std$kcdv-mean(POLICYII.std$kcdv,na.rm=T))/sd(POLICYII.std$kcdv,na.rm=T)
POLICYII.std$nopolicy<-(POLICYII.std$nopolicy-mean(POLICYII.std$nopolicy,na.rm=T))/sd(POLICYII.std$nopolicy,na.rm=T)

###STANDARDIZED RANDOM EFFECTS MODELS OF THE DISASTER'S EFFECT ON AVERAGE###
std1=plm(dependent~psdv-1,data=POLICYII.std,index=c("city","year"),model="random");summary(std1)

std2=plm(dependent~psdv*npo-1,data=POLICYII.std,index=c("city","year"),model="random");summary(std2)

std3=plm(dependent~psdv*enpo-1,data=POLICYII.std,index=c("city","year"),model="random");summary(std3)

std4=plm(dependent~psdv*gnpo-1,data=POLICYII.std,index=c("city","year"),model="random");summary(std4)

###STANDARDIZED RANDOM EFFECTS MODELS OF THE DISASTER'S EFFECT CONTINGENT ON POLICY###
std5=plm(dependent~ttdv+kcdv+nopolicy-1,data=POLICYII.std,index=c("city","year"),model="random");summary(std5)

std6=plm(dependent~ttdv*npo+kcdv*npo+nopolicy*npo-1,data=POLICYII.std,index=c("city","year"),model="random");summary(std6)

std7=plm(dependent~ttdv*enpo+kcdv*enpo+nopolicy*enpo-1,data=POLICYII.std,index=c("city","year"),model="random");summary(std7)

std8=plm(dependent~ttdv*gnpo+kcdv*gnpo+nopolicy*gnpo-1,data=POLICYII.std,index=c("city","year"),model="random");summary(std8)

###POST HOC T-RATIOS###
#Model 2
#Minimium
effect<-std2$coefficients[1]+(std2$coefficients[3]*min(POLICYII.std$npo)); effect
var<-std2$vcov[1,1]+2*min(POLICYII.std$npo)*std2$vcov[1,3]+(min(POLICYII.std$npo)^2)*std2$vcov[3,3]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-std2$coefficients[1]+(std2$coefficients[3]*max(POLICYII.std$npo)); effect
var<-std2$vcov[1,1]+2*max(POLICYII.std$npo)*std2$vcov[1,3]+(max(POLICYII.std$npo)^2)*std2$vcov[3,3]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 3
#Minimium
effect<-std3$coefficients[1]+(std3$coefficients[3]*min(POLICYII.std$enpo)); effect
var<-std3$vcov[1,1]+2*min(POLICYII.std$enpo)*std3$vcov[1,3]+(min(POLICYII.std$enpo)^2)*std3$vcov[3,3]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-std3$coefficients[1]+(std3$coefficients[3]*max(POLICYII.std$enpo)); effect
var<-std3$vcov[1,1]+2*max(POLICYII.std$enpo)*std3$vcov[1,3]+(max(POLICYII.std$enpo)^2)*std3$vcov[3,3]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 4
#Minimium
effect<-std4$coefficients[1]+(std4$coefficients[3]*min(POLICYII.std$gnpo)); effect
var<-std4$vcov[1,1]+2*min(POLICYII.std$gnpo)*std4$vcov[1,3]+(min(POLICYII.std$gnpo)^2)*std4$vcov[3,3]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-std4$coefficients[1]+(std4$coefficients[3]*max(POLICYII.std$gnpo)); effect
var<-std4$vcov[1,1]+2*max(POLICYII.std$gnpo)*std4$vcov[1,3]+(max(POLICYII.std$gnpo)^2)*std4$vcov[3,3]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 6
#Minimium
effect<-std6$coefficients[1]+(std6$coefficients[5]*min(POLICYII.std$npo)); effect
var<-std6$vcov[1,1]+2*min(POLICYII.std$npo)*std6$vcov[1,5]+(min(POLICYII.std$npo)^2)*std6$vcov[5,5]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-std6$coefficients[1]+(std6$coefficients[5]*max(POLICYII.std$npo)); effect
var<-std6$vcov[1,1]+2*max(POLICYII.std$npo)*std6$vcov[1,5]+(max(POLICYII.std$npo)^2)*std6$vcov[5,5]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 7
#Minimium
effect<-std7$coefficients[1]+(std7$coefficients[5]*min(POLICYII.std$enpo)); effect
var<-std7$vcov[1,1]+2*min(POLICYII.std$enpo)*std7$vcov[1,5]+(min(POLICYII.std$enpo)^2)*std7$vcov[5,5]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-std7$coefficients[1]+(std7$coefficients[5]*max(POLICYII.std$enpo)); effect
var<-std7$vcov[1,1]+2*max(POLICYII.std$enpo)*std7$vcov[1,5]+(max(POLICYII.std$enpo)^2)*std7$vcov[5,5]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

#Model 8
#Minimium
effect<-std8$coefficients[1]+(std8$coefficients[5]*min(POLICYII.std$gnpo)); effect
var<-std8$vcov[1,1]+2*min(POLICYII.std$gnpo)*std8$vcov[1,5]+(min(POLICYII.std$gnpo)^2)*std8$vcov[5,5]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p
#Maximum
effect<-std8$coefficients[1]+(std8$coefficients[5]*max(POLICYII.std$gnpo)); effect
var<-std8$vcov[1,1]+2*max(POLICYII.std$gnpo)*std8$vcov[1,5]+(max(POLICYII.std$gnpo)^2)*std8$vcov[5,5]; var
se<-sqrt(var); se
t<-effect/se; t
p<-(1-pt(abs(t),102))*2;p

